import plotly
import plotly.graph_objs as go
import plotly.express as px
from plotly.subplots import make_subplots
import torch
import numpy as np
import pandas as pd
import dill as pickle
import seaborn as sns
import random
import math
Wave equation with one spatial variable
$$ \frac{\partial^{2} u}{\partial t^{2}} - \frac{1}{25} \frac{\partial^{2} u}{\partial x^{2}} = 0, $$ $$ \\ 71\times71, x \in [0; 1], t \in [0; 1]. $$
grid_res = 70
title = 'wave_equation'
df = pd.read_csv('C:\\Users\\Ksenia\\NSS\\ODE_projects\\SDE-DS\\wave_sln_70.csv', header=None)
initial_data = torch.from_numpy(df.to_numpy()).reshape(-1)
x_grid = np.linspace(0, 1, grid_res + 1)
t_grid = np.linspace(0, 1, grid_res + 1)
params = [x_grid, t_grid]
x = torch.from_numpy(x_grid)
t = torch.from_numpy(t_grid)
grid = torch.cartesian_prod(x, t).float()
u_sde = torch.load('C:\\Users\\Ksenia\\NSS\\ODE_projects\\SDE-DS\\average_wave_pt\\sols_stacked_35_solutions_shape_torch.Size([5041, 1]).pt', weights_only=True)
u_sde_mean_tens = u_sde.reshape(*[len(i) for i in params]).detach().cpu().numpy().T
u_main = torch.load('C:\\Users\\Ksenia\\NSS\\ODE_projects\\SDE-DS\\file_u_main_[35, 71, 71]_autograd_[0].pt', weights_only=False)
u_main.shape
(35, 71, 71)
def calculate_statistics(u_main):
mean_arr = np.zeros((u_main.shape[1], u_main.shape[2]))
var_arr = np.zeros((u_main.shape[1], u_main.shape[2]))
s_g_arr = np.zeros((u_main.shape[1], u_main.shape[2])) # population standard deviation of data.
s_arr = np.zeros((u_main.shape[1], u_main.shape[2])) # sample standard deviation of data
for i in range(u_main.shape[1]):
for j in range(u_main.shape[2]):
mean_arr[i, j] = np.mean(u_main[:, i, j])
var_arr[i, j] = np.var(u_main[:, i, j])
m = np.mean(u_main[:, i, j])
s_arr[i, j] = math.sqrt(np.sum(list(map(lambda x: (x - m)**2, u_main[:, i, j])))/(len(u_main[:, i, j]) - 1))
mean_tens = torch.from_numpy(mean_arr)
var_tens = torch.from_numpy(var_arr)
s_g_arr = torch.from_numpy(var_arr) ** (1/2)
s_arr = torch.from_numpy(s_arr)
# Confidence region for the mean
upper_bound = mean_tens + 1.96 * s_arr / math.sqrt(len(u_main))
lower_bound = mean_tens - 1.96 * s_arr / math.sqrt(len(u_main))
mean_tens = mean_tens.reshape(-1)
upper_bound = upper_bound.reshape(-1)
lower_bound = lower_bound.reshape(-1)
return mean_tens, upper_bound, lower_bound
u_bs_mean_tens, u_bs_upper_bound, u_bs_lower_bound = calculate_statistics(u_main)
Метод Соболя основан на разложении дисперсии функции $f(\mathbf{x})$ на вклады различных входных параметров. Формула выглядит следующим образом:
$$ f(\mathbf{x}) = f_0 + \sum_{i=1}^n f_i(x_i) + \sum_{1 \leq i < j \leq n} f_{ij}(x_i, x_j) + \cdots + f_{1,2,\ldots,n}(x_1, x_2, \ldots, x_n), $$ где:
SDE: Представляет компоненты, связанные с дисперсией первого порядка, т.е. с влиянием индивидуальных стохастических параметров.
BS: Преимущественно соответствует дисперсиям второго порядка, отражая условные вероятности и зависимости между параметрами.
$$ M[BS−M[SDE]] $$ BS — результат применения подхода на основе байесовской сети.
SDE — результат применения подхода на основе стохастических дифференциальных уравнений.
u_temp = u_main - u_sde_mean_tens
result_1 = np.mean(u_temp , axis=0)
print(result_1)
result_1.shape
[[-0.01292047 -0.01334606 -0.01373914 ... 0.01158872 0.0090557 0.0050815 ] [ 0.04311887 0.0433047 0.04350033 ... -0.16430853 -0.1725123 -0.18118833] [ 0.08509692 0.085596 0.08606815 ... -0.29792038 -0.31092852 -0.32398838] ... [ 0.14299496 0.14392528 0.1447486 ... -0.5721952 -0.5990019 -0.6262118 ] [ 0.07228713 0.07265554 0.07298737 ... -0.29356402 -0.31186715 -0.3312171 ] [-0.01658184 -0.01707336 -0.01756203 ... 0.03265116 0.0267742 0.01848607]]
(71, 71)
# building 3-dimensional graph
fig = go.Figure(data=[
go.Mesh3d(x=grid[:, 0], y=grid[:, 1], z=initial_data, name='Initial field',
legendgroup='i', showlegend=True, color='rgb(139,224,164)',
opacity=0.5),
go.Mesh3d(x=grid[:, 0], y=grid[:, 1], z=result_1.reshape(-1), name='M[BS−M[SDE]]',
legendgroup='s', showlegend=True, color='lightpink',
opacity=1),
go.Mesh3d(x=grid[:, 1], y=grid[:, 1], z=u_sde_mean_tens.reshape(-1), name='Solution field - M[SDE]',
legendgroup='sde', showlegend=True, color='rgb(139,200,164)',),
])
fig.update_layout(scene_aspectmode='auto')
fig.update_layout(showlegend=True,
scene=dict(
xaxis_title='x1',
yaxis_title='x2',
zaxis_title='u',
zaxis=dict(nticks=10, dtick=1),
aspectratio={"x": 1, "y": 1, "z": 1}
),
height=800, width=800
)
fig.show()
$$ M[BS]−M[SDE] $$
u_main_bs = np.mean(u_main, axis=0)
result_2 = u_main_bs - u_sde_mean_tens
print(result_2)
result_2.shape
[[-0.01292047 -0.01334606 -0.01373914 ... 0.01158872 0.0090557 0.00508148] [ 0.04311887 0.04330469 0.04350033 ... -0.16430855 -0.17251235 -0.18118832] [ 0.08509695 0.08559598 0.08606817 ... -0.29792035 -0.31092864 -0.32398844] ... [ 0.14299497 0.1439253 0.14474854 ... -0.57219523 -0.5990019 -0.62621176] [ 0.07228713 0.07265553 0.07298738 ... -0.29356405 -0.31186718 -0.33121714] [-0.01658184 -0.01707336 -0.01756203 ... 0.03265116 0.02677419 0.01848608]]
(71, 71)
rmse = np.sqrt(np.mean((result_1 - result_2) ** 2))
print(f"Root Mean Squared Error (RMSE): {rmse}")
Root Mean Squared Error (RMSE): 3.5406208098720526e-07
# building 3-dimensional graph
fig = go.Figure(data=[
go.Mesh3d(x=grid[:, 0], y=grid[:, 1], z=initial_data, name='Initial field',
legendgroup='i', showlegend=True, color='rgb(139,224,164)',
opacity=0.5),
go.Mesh3d(x=grid[:, 0], y=grid[:, 1], z=result_1.reshape(-1), name='M[BS−M[SDE]]',
legendgroup='s', showlegend=True, color='lightpink',
opacity=1),
go.Mesh3d(x=grid[:, 0], y=grid[:, 1], z=result_2.reshape(-1), name='M[BS]−M[SDE]',
legendgroup='c', showlegend=True, color='blue',
opacity=1),
go.Mesh3d(x=grid[:, 1], y=grid[:, 1], z=u_sde_mean_tens.reshape(-1), name='Solution field - M[SDE]',
legendgroup='sde', showlegend=True, color='rgb(139,200,164)',),
])
fig.update_layout(scene_aspectmode='auto')
fig.update_layout(showlegend=True,
scene=dict(
xaxis_title='x1',
yaxis_title='x2',
zaxis_title='u',
zaxis=dict(nticks=10, dtick=1),
aspectratio={"x": 1, "y": 1, "z": 1}
),
height=800, width=800
)
fig.show()
max_value = np.max(result_1)
variance = np.var(u_sde_mean_tens)
print(f"max M[BS - M[SDE]]: {max_value}")
print(f"var SDE: {variance}")
max M[BS - M[SDE]]: 0.4558561146259308 var SDE: 3.770282030105591
x = np.linspace(-5, 5, 100)
y = np.linspace(-5, 5, 100)
fig = make_subplots(rows=1, cols=2, subplot_titles=('<b>Solution<b>', '<b>M[SDE]<b>'))
fig.add_trace(
go.Heatmap(
z=u_sde_mean_tens, x=x, y=y, colorscale='magma', opacity=1,
showscale=True,
colorbar=dict(x=1.0)
),
row=1, col=1
)
fig.add_trace(
go.Heatmap(
z=u_sde_mean_tens, x=x, y=y, colorscale='magma', opacity=1,
showscale=False
),
row=1, col=2
)
font=dict(
family="ATimes New Roman, Times, serif",
size=16,
color="black"
),
title=dict(
font=dict(
family="Times New Roman, Times, serif",
size=20,
color="black"
)
)
fig.show()
fig = make_subplots(
rows=1, cols=1,
subplot_titles=('<b>M[BN]−M[SDE]</b>'))
fig = go.Figure(
data=go.Heatmap(
z=result_2,
x=x,
y=y,
colorscale='magma',
opacity=1,
showscale=True,
colorbar=dict(x=1.02)
)
)
fig.update_layout(
title={
'text': "<b>M[BN]−M[SDE]</b>",
'y': 0.9,
'x': 0.5,
'xanchor': 'center',
'yanchor': 'top', },
height=600,
width=600,
showlegend=False,
xaxis=dict(scaleanchor="y"),
yaxis=dict(scaleanchor="x"),
xaxis2=dict(scaleanchor="y2"),
yaxis2=dict(scaleanchor="x2")
)
font=dict(
family="ATimes New Roman, Times, serif",
size=16,
color="black"
),
title=dict(
font=dict(
family="Times New Roman, Times, serif",
size=20,
color="black"
)
)
fig.show()